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ABSTRACT 

Context. Massive stars play a vital role in the Universe. However, their evolution even on the main sequence is not yet well understood. 
Aims. Due to the steep mass-luminosity relation, massive main sequence stars become extremely luminous. This brings their envelopes 
very close to the Eddington limit. We are analysing stellar evolutionary models in which the Eddington limit is reached and exceeded, 
and explore the rich diversity of physical phenomena which take place in their envelopes, and we investigate their observational 
consequences. 

Methods. We use the published grids of detailed stellar models by Brott et al. (2011) and Kohler et al. (2015), computed with a 
state-of-the-art one-dimensional hydrodynamic stellar evolution code using EMC composition, to investigate the envelope properties 
of core hydrogen burning massive stars. 

Results. We find that at the stellar surface, the Eddington limit is almost never reached, even for stars up to 500 Mq. When we define 
an appropriate Eddington limit locally in the stellar envelope, we can show that most stars more massive than ~ 40 Mq actually exceed 
this limit, in particular in the partial ionization zones of iron, helium or hydrogen. While most models adjust their structure such that 
the local Eddington limit is exceeded at most by a few per cent, our most extreme models do so by a factor of more than seven. We 
find that the local violation of the Eddington limit has severe consequences for the envelope structure, as it leads to envelope inflation, 
convection, density inversions and possibly to pulsations. We find that all models with luminosities higher than 4 x 10^ Lq, i.e. stars 
above ~ 40 Mq show inflation, with a radius increase of up to a factor of about 40. We find that the hot edge of the S Dor variability 
region coincides with a line beyond which our models are inflated by more than a factor of two, indicating a possible connection 
between S Dor variability and inflation. Furthermore, our coolest models show highly inflated envelopes with masses of up to several 
solar masses, and appear to be candidates to produce major LBV eruptions. 

Conclusions. Our models show that the Eddington limit is expected to be reached in all stars above ~ 40 Mq in the LMC, and 
by even lower mass stars in the Galaxy, or in close binaries or rapid rotators. While our results do not support the idea of a direct 
super-Eddington wind driven by continuum photons, the consequences of the Eddington limit in the form of inflation, pulsations and 
possibly eruptions may well give rise to a significant enhancement of the the time averaged mass loss rate. 
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1. Introduction 


Massive stars are powerful engines and strongly affect the 
evolution of star fo rming galaxies throughout cosmic time 
dBresolin et al.|[200^ . In particular the mos t massive ones pro - 
duce copious amounts of ionisin g photons (jPoran et^ [201^ 
emit p owerful stellar winds (e.g.. ]Kudritzki & Pulsll200C : ISmithI 
l2014li and in their final explosions are suspected to produce 
the most energetic an d spectacular stellar explosions, as pair- 
instability supemoyaej^z^revaeTaOlQT^ij superluminous su- 
pernovaejGaL^rmet^ 2009h. and long-dur ation gamma-ray 
bursts (iLarsson et al.ll200'^ Raskin et al.ll2008l) . 

Massive main sequence stars, which we understand here as 
those which undergo core hydrogen burning, have a much higher 
luminosity than the Sun, as they are known to obey a sim¬ 
ple mass-luminosity relation, L ~ M", with a > 1. However, 
whereas thi s relation is very steep near the Solar mass (a ^ 5), it 
is shown in Kippenha hn & W eigerll (1 19901) that cr —> 1 for M 
oo. Indeed. iKdhler et al.l (|2015l) find a ^ 1.1 forM = 500 Mq. 

Since the Eddington factor is proportional to LjM, it is de¬ 
bated in the literature whether main sequence stars of hig her and 
higher initial mass eventually reach the Eddington-limit (iLangej 
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ICrowther et alll2010l: iMaederet al.ll2012l) . The answer is 
clearly; yes, they do. Even when only electron scattering is con¬ 
sidered as a source of radiative opacity, the Eddington-limit cor¬ 
responds to a luminosity-to-mass ratio of 7? ;= log(^/^) ^ 
4.6 (iLanger & Kudritz^l2014l) for hot stars with a solar helium 
abundance. This is extremely close to the !??-values obtained for 


models of su permassive stars, where this ratio is nearly mass- 
independent (lEuller et alJl986l:l^oll986h . Infact. lKatol(ll986l) 
showed that zero-age main sequence models computed only with 
electron scattering opacity do reach the Eddington limit at a mass 
of about ~ 10^ Mq. 

Whether supermassive stars exist is an open question. 
Also the mass limit of ordinary stars is presently uncertain 
(ISchneider et al.l l2Q14t) . However, there is ample evidence for 
stars with initial masses well above 100 Mq in the local Universe. 
A number of close binary stars hav e been found with compo¬ 
nent initial masses above 100 M^ (Schnurr et al.l 200^ 120091 : 
iTavlor et^l201 UlSana et al.ll2013l) . Crowther et al.l ( ~1 Ol) pro- 
posed initial masses of up to 300 Mq for several stars in the 
Large Magellanic Cloud (LMC), based on their luminosities. 
iBestenlehner et all (l2014h identified more than a dozen stars 
more massive than IOOMq from the sample of ~ 1000 OB stars 
near 30 Doradus, which are analysed i n the frame of VLT-Elames 
Tarantula Survey (lEvans et alJl201 ih . The hydrogen-rich stars 
among them have measured 7?-values of up to 4.3. In hot stars 
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Fig. 1. The different Eddington factors inside a 285 Mq, non¬ 
rotating, main sequence model with logL/Lo = 6.8 and 
Teff = 46600 K (cf. Fig. |9] and Appendix 1^. The shaded areas 
mark the different convection zones and the hatched area marks 
the region with a density inversion. The radius of the un-inflated 
core is denoted as (defined in Sect. S. The black dashed 
horizontal line is drawn at T = 1 for convenience. 


with hnite metallicity, the i on opacities can easily exceed the 
electron scattering opacity (llglesias & Rogerslll996l) . It is thus 
to be expected that the true Eddington limit, which accounts for 
all opacity sources, is located at '??-values of 4.3 or below. This 
implies that these stars should in fact also have reached, or ex¬ 
ceeded their Eddington limit. 

In this paper, we explore the question of massive main se¬ 
quence stars reaching, or exceeding the Eddington limit from 
the theoretical side. We show by means of detailed stellar mod¬ 
els as described in Sectionl^that even all stars more massive than 
~ 40 Mq are found to reach the Eddington limit. In Section[3l we 
demonstrate the need to properly define a local Eddington factor 
in the stellar interior, which we then use in Section]?] to show 
that when it exceeds the critical value of one, the stellar enve¬ 
lope becomes inflated. We show further in Section]5] and]^ that 
super-Eddington conditions can lead to density inversions, and 
induce convection. We compare our results to previous studies 
in Section]?) and relate them to observations in Section]8] before 
summarising our conclusions in Section]9] 

2. Stellar models 

The grids of stellar models use d for the present s t udy ha ve been 
published in iBrott et al.l (1201 ll) and ]K6hler et alJ (l2()15l) . In this 
paper, we consider only the core hydrogen burning models com¬ 
puted with E MC metallicity. E ach s tellar evolution seq uence 
computed by iBrott et ^ (1201 Ih and iKohler et al.l (1201 5l) con¬ 
sists of typically 2000 individual stellar models. However, the 
full amount of data dehning a stellar model is only stored for a 
few dozen time points per sequence, in non-regular intervals. It 
is those stored models which are analysed here. This scheme has 
the disadvantage that the density of models in the investigated 
parameter space is not always as high as it should be ideally. 
Still, as shown below, it allows for a thorough sampling of the 
considered parameter space, and it is fully consistent with the 
results already published. 

The stellar models were computed with a state-of-the-art 
one-dimensional hydrodynamic implicit Lagrangian code (BEC) 


which incorporates latest input physics (for details, see Braunl 

ll997l:lYoon et alJ2006l:]Brott et al.l201 lllKohler et alJ2015L and 

references therein). Convection was trea ted as in the standard 
non-adiabatic mixing length approach (iBohm-Vitensd 119581 : 
iKippenhahn & W eiger 3]T^ a nd a mixing length parameter of 
a - ////p = 1.5 (ILangei 1199 Ih was adopted, with I and Hp be¬ 
ing the mixing length and the pressure scale height respectively. 
This value of the mixing l ength parameter d oes lead to a good 
representation of the Sun (ISuiis et alJl2008l) . whereas its cali¬ 
bration to multi-dimensional hydrodynamic models shows that 
it tends to decrease towa rds lower gravities (iTramnedach et al.l 
120141 : ]Magic et all 1201 5ll . The convective velocities were lim¬ 
ited to the local value of the a diabatic sound speed. The con¬ 
tribution of turbulent pressure (Ide Jagei]ll984l) was neglected, 
since it is not expecte d to be importan t in determining the stellar 
hydrostatic structure (iStothersll2003l) . Indeed, our recent study 
which includes turbulent pressure (Grassitelli et al. in prepara¬ 
tion) shows that e.g. for an 80 Mq evolutionary sequence, the 
stellar radius is increased over that of models without turbu¬ 
lent pressure by at most a few per cents at any time during 
its main-sequence evolution. Rotation al mixing of chemical el¬ 
ements following iHeger et"^ (i2000h and transport of angular 
momentum by mag netic fields du e to the Spruit-Taylor dynamo 
were also included (ISpruitll2002h . The efficiency parameters fc 
and fp for rotational mix ing were set to 0.0228 and 0.1 respec¬ 
tively (iBrott et ^1201 1]). Radiative opacities w ere interpolated 
from the OPAL tables (llglesias & Rogersll996l) . The opacity en¬ 
hancement due to Fe-group elements at T ~ 200 kK plays a vital 
role in determining the envelope structure in our stellar models. 
We note that even though flux-mean opacities are appropriate 
to study the momentum balance near the stellar photosphere, 
we only consider the Rosseland mean opacities in the follow¬ 
ing, which are thought to behave very similarly to the flux-mean 
opacities especially at an optical depth larger than one. 

The outer boundary condition of the stellar models corre¬ 
sponds to a plane-parallel gray atmosphere model on top of the 
photosphere. In other words, the effective temperature was used 
as a boundary condition at a Rosseland optical depth of 2/3. The 
adopted stellar wind mass loss recipe does lead to small but fi¬ 
nite outflow velocities in the outermost layers, which induces a 
slight deviation from hydrostatic equil ibrium. 

The mass loss prescription from IVink et all (l2000i 1200 ih 
was employed to account for the winds of O- and B- 
fype stars. Moreover, parameter ized mass loss rates from 
iNieuwenhuiizen & de .lageil (Il990h were used on the cooler side 
of the bi-stabilit y jump, i.e. at effective temperature s less than 
22000 K, if the iNie uwenhuiizen & de Jagei) (Il990ll mass loss 
rate exceeded that of l Vink et alfil2000l]200llk Wolf-Rayet (WR) 
type mass los s was accounted for using the empirical pre- 
scription fro m Hamann et alJ (1 19951) divided by a factor of 10 
(I Yoon et ^120061) . when the surface helium mass fraction be¬ 
came greater than 70%. 

Evolutionary sequences of massive stars, with and without 
rotation, were computed up to an initial mass of 500 Mq, starting 
with EMC composition. The initial mass fractions of hydrogen, 
helium and metals were taken to be 0.7391, 0.2562 and 0.0047 
respectively, in accord ance with the obs ervations of young mas¬ 
sive stars in the EMC (iBrott et alJl201 ih . 


3. The Eddington Limit 

The Eddington limit refers to the condition where the outward 
radiative acceleration in a star balances the inward gravity, in 
hydrostatic equilibrium. It is a concept which is thought to ap- 
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Fig. 2. Positions of the analysed stellar models with Fmax > 0.9 in the Hertzsprung-Russell diagram (colore d dots). Models wit h 
Tmax > 1-1 are colored dark-blue. The solid lines show the evolutionary tracks of non-rotating stellar models (iKohler et al.ll2015h . 
The initial masses are marked in units of solar mass. The dashed line c orresponds to the z ero-age main sequence of the non-rotating 
models. The hot and the cool edges of the S Dor instability strip from (ISmith et al.ll2004h are indicated with thick dotted lines. The 
interior structures of the models marked with yellow diamonds are shown in AppendixlDl 


ply at the stellar surface, in the sense that if the Eddington 
limit is exceeded, a mass outflow should arise dEddingtonl 1 92(i 
lOwocki et al.ll2004l) . If we denote the gravity as g - GMjr^ and 
the radiative acceleration mediated through the electron scatter¬ 
ing opacity as then the classical Eddington fac¬ 

tor Eg is defined as 

HoJ—i 

Eg ^ (1) 

Tfidd g A-ttcGM 

where L,M and are the luminosity, mass and electron¬ 
scattering opacity respectively, with the physical constants hav¬ 
ing their usual meaning. The classical Eddington parameter Eg 
therefore does not depend on the radius r as the inverse scal¬ 
ing in both 0rad and g cancel out. Whereas Eg is often convenient 
to consider, it provides a sufficient instability criterion to stars 
but not a necessary one, because usually the true opacity does 
exceed the the electron scattering opacity significantly and also 
contributes to the radiative force. 

As it turns out below, even when the Rosseland mean opaci¬ 
ties are used, the models analysed in this paper practically never 
reach the Eddington limit at their surface. Therefore, we instead 


consider the Eddington factor in the stellar interior as 

ETrl ^ 

l-Eddir) 4ncGM(ry 


( 2 ) 


where M(r) is the Lagrangian mass coordinate, K(r) is the 
Rosse land mean opacity and Ur) is the local luminosity (iLanged 
Il997h . However, r'(r) > 1 also does not provide a stability limit 
in the stellar interior because the stellar layers turn convectively 
unstable following Schwarzsc hild’s criterion when r'(r) ^ 1 
(iJoss et al.ll973tlLangedl9^ . As the luminosity transported by 
convection does not contribute to the radiative force, we subtract 
the convective luminosity in the above expression and redefine 
the Eddington factor as 


TradCO _ TcQnv(r)) 

Uddir) AncGM{r) 


(3) 


Eor example, near the stellar core where convective energy trans¬ 
port is highly efficient, E(r) stays well below unity in spite of 
E'(r) » 1 and no instability, i.e. departure from hydrostatic equi¬ 
librium, occurs (see Eig. [1]). In the rest of the paper we will refer 
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to r(r) as the Eddington factor unless explicitly specified other¬ 
wise. 

Even with this dehnition, a super-Eddington layer inside a 
star does not necessarily lead to a departure from hydrostatic 
equilibrium or a sustained mass outflow. In the outer envelopes 
of massive stars non-adiabatic conditions prevail and convective 
energy transport is highly inefficient which pushes r(r) close 
to (or above) one. We find that the stellar models counteract 
such a super-Eddington luminosity by developing a posi tive gas 
pressure gradient, thu s restoring hydrostatic equilibrium (iLangej 
ll997HAsr)lund|[T99^ . In such situations the canonical definition 
of Ledd being the maximum sustainable radiative luminosity lo¬ 
cally in the stellar interior (in hydrostatic equilibrium) breaks 
down and loses its significance. As we shall see below the radia¬ 
tive luminosity beneath the photosphere can be up to a few times 
the Eddington luminosity. 

In Eig. [1] the behavior of E and E' is shown along with 
the electron-scattering Eddington factor Ee in a 285 Mq non¬ 
rotating stellar model, which provides an educative example (see 
Appendix iP] for further examples). As explained above, E' and 
Ee are significantly greater than one in the convective core of the 
star. The indicated sub-surface convection zones are caused by 
the opacity peaks at E ~ 1.5 x 10® K (deep iron bump) and at 
T ~ 2 X 10® K (iron bump). Near the bottom of the inflated 
envelope {rfrcoK ^ 1; see Sect. |4] for the definition of rcore), T 
approaches one and the Ee opacity bump drives convection. An 
extended region with E 1 follows. A thin shell very close to the 
photosphere contains the layers with a positive density gradient 
and with E > 1. 



Fig. 4. Maximum Eddington factor (Emax, red line), effective 
temperature (Eefr, blue line) and effective temperature at the non- 
inflated core (Eefr^core, pink line, see Sect.|4] ean.fTOl) as a function 
of initial mass for the non-rotating ZAMS models. The black 
dotted line at Emax = 1 is drawn for reference. 

namic time scale reveal that some, and potentially many, of our 
models are pulsationally unstable, as will be shown in a forth¬ 
coming paper. However, in the cases analysed so far, the pulsa¬ 
tions saturate and do not lead to a destruction or ejection of the 
inflated envelopes. In this respect, we consider the analysis of 
the hydrostatic equilibrium structures performed in this paper as 
usefuL 


E 

U 


60 50 40 30 20 10 0 

Teff [kK] 

Fig. 3. Maximum Eddington factor in the analysed stellar mod¬ 
els as a function of Teff for all stellar models shown in Eig. |2] 
i.e. with Emax > 0.9. The three peaks in Emax at Teff/kK of ~ 
55, 25 and 5.5 correspond to the three opacity bumps associated 
with the ionization zones of Ee, Hell and H respectively. Inset. 
Models with E^ax > 0.94 and Teff/kK > 10 kK. The blue hori¬ 
zontal line at E^ax = 1 is drawn for reference. 



We note that the stellar models have been computed with 
a hydrodynamic stellar evolution code. However, due to the 
large time steps required for stellar evolution calculations, non¬ 
hydrostatic solutions are suppressed by our numerical scheme. 
The resulting hydrostatic structures are still val i d solu tions 
of the hydrodynamic equations (see iHeger et^ (l200(]l) and 
iKozvreva et al.l (12014h for the equations to be solved). Models 
computed with time steps small enough to resolve the hydrody- 


3.1. Effect of rotation on the Eddington iimit 


The effect of the centrifugal force on the structure of rotat¬ 
ing stellar models has been stu d ied by a n umber of groups in 
the pa st, including iHeger et"^ (l2000l) and iMaeder & Mevnej 
(l2000h . This is done by describing the models in a 1-D ap¬ 
proximation w here all variables are taken as averages over iso- 
baric surfaces (iKippenhahn & Thomasl[T970ll . The stellar struc¬ 
ture equation s are modified to in clude the effect of the cen¬ 
trifugal force (lEndal & Sofial[T97^ . The equation of hydrostatic 
equilibrium becomes 


dP 

dm 


47Tr^ + fp 


GM{r) 


= 0 , 


(4) 


and the radiative temperature gradient in the energy transport 
equation (in the absence of convection) takes the form 


Vrad - 


3 kPL fr 
IhnacG MT* fp’ 


(5) 


w here the quantities fp and fp have the same definition as 
in [Heger et al.l (120001) . Consequently, the Eddington luminosity 
gets modified as: 


Tudd - 


AjtcGM fp 
K ft 


(6) 


However, the Eddington factor. 


Li-ad _ _V_^ _ Aa E^V 

Tudd ^rad 3 P 


does not have any explicit dependence on fp and fp because the 
factor fpj fp cancels out. Therefore formally, the Eddington fac¬ 
tor remains unaffected by rotation. Of course, if the internal evo¬ 
lution of a rotating model is changed, for example by rotational 


4 















































D. Sanyal et al.: Massive main sequence stars evolving at the Eddington limit 



60 50 40 30 20 10 0 


Teff [kK] 


Fig. 5. Analyzed models shown as colored dots on the spectroscopic HR diagram (sHRD) with the color representing the value of 
Fmax in each model. The Eddington factor Ee (assuming electron scattering opacity only) is directly proportional to the quantity 
jSf and is indicated on the right Y-axis for a hydrogen mass fraction of Y = 0.74. Some representative evolutionary tracks of non¬ 
rotating models, for different initial masses (indicated along the tracks in units of solar mass), are also shown with solid black lines. 
The green and blue dotted horizontal lines correspond to Eg = 0.7 for X = 0.7 and Y = 0 respectively. The color palette and the 
models marked with yellow diamonds correspond to those in Eig.|3 The black dashed line roughly divides the sHRD into distinct 
regions with Tmax > 1 and Tmax < 1 ■ 


mixing, its Eddington factor will still be different from that of 
the corresponding non-rotating model. 

Of course, real stars are three-dimensional and the centrifu¬ 
gal force must affect the hydrostatic stability limit. However, this 
is expected to be a function of the latitude at the stellar surface, 
and in a 2-D view, the effect will be largest at the equator dLangej 
Il997h . To hrst order, the critical luminosity Lg to unbind matter 
at the stellar equatorial surface becomes 


Eg — TjJdd 


ttrot 

^Kep 


2 \ 


(8) 


where Uiot and UKgp are the stellar equatorial rotation velocity 
and the corresponding Keplerian value, respectively. However, 
to compute the effect reliably, the stellar deformation due to ro¬ 
tation as we ll as the effect of gravity darkening nee d to be ac¬ 
counted for dMaeder & MevnefcoOOtlMaederllTOO^ . To do this 
realistically for stars near the Eddington limit requires at least 
2-D calculations. 

The implication is that the effect of rotation on the criti¬ 
cal stellar luminosity can not be properly described through the 


models analysed here. Those models see the same critical lumi¬ 
nosity as if rotation was absent. Since mixing of helium in these 
models is very weak for rotation rates below the ones required 
for chemically homogeneous evolution, most of t he rotating 
models evolve very simil ar to the non-rotating ones dBrott et al.l 
1201 ibf^hler et al]l20L5h . and thus merely serve to augment our 
database. 


3.2. The maximum Eddington factor 

In our stellar models, we have determined the maximum 
Eddington factor Emax over the whole star, i.e Emax 
max [E(r)]. The maximum Eddington factor Emax generally oc¬ 
curs in the outer envelopes of our models, where convective en¬ 
ergy transport is much less efficient than in the deep interior. 
The variation of E^ax across the upper HR diagram is shown in 
Eig.|2]for all analysed core hydrogen burning models which have 
Fmax > 0.9. 

Three distinct regions with E^ax > 1 can be identihed in 
Eig.|2] which can be connected to the opacity peaks of iron, he¬ 
lium and hydrogen. When one of these opacity peaks is situated 
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TefflkK] 

Fig. 6. Hertzspmng-Russell diagram showing the logarithm of 
the optical depth t at the position of Fmax in color, for all anal¬ 
ysed models which have Fmax > 0.9. Some representative evolu¬ 
tionary tracks of non-rotating models, for different initial masses 
(indicated along the tracks in units of solar mass), are also shown 
with solid black lines. 



60 50 40 30 20 10 0 

T,i[kK] 

Fig. 7. The Eddington factors at the stellar surface, F(R*), are 
shown on the sHRD for our analysed set of models. The models 
with F(R*) < 0.9 are indicated with smaller black circles and 
those with F(R*) > 0.9 are shown as colored dots. 


sufficiently close to the stellar photosphere, the densities in these 
layers are so small that convective energy transport becomes 
inefficient. As a consequence, super-Eddington layers develop 
which are stabilized by a positive (i.e. inward directed) gradient 
in density and gas pressure (see Sect.|5]below). The envelope in¬ 
flation which occurs when Fmax approaches one is discussed in 
Sect.|4] 

Eigure |3] shows Fmax as a function of the effective tempera¬ 
ture for all our models with Fmax > 0.9. The models which have 
the hydrogen opacity bump close to their photosphere can obtain 
values of Fmax as high as ~ 7. This manifests itself as a promi¬ 
nent peak around Teff ^5.5 kK. The inset panel shows the much 
weaker peaks in Fmax due to the partial ionization zones of Fe 
and Hell, at T/kK ~ 200 and 50 respectively. The peak caused 
by the Fe opacity bump may extend to hotter effective temper¬ 
atures and apply to hot, hydrogen-free Wolf-Rayet stars, which 
are not part of our model grid. 


For stars above about 125 Mq, Fmax reaches one, even on the 
zero-age main sequence. This is demonstrated in Fig. 0 ] which 
shows both Fmax and effective temperature as a function of mass 
for the non-rotating stellar models. As these models evolve away 
from the ZAMS to cooler temperatures, super-Eddington layers 
develop in their interior. The blue curve shows a maximum ef¬ 
fective temperature of 57 000 K at about 200 Mq, beyond which 
it starts decreasing with further increase in mass. This behaviour 
is related to the phenomenon of inflation which is discussed in 
detail in Sect.|4] However, the ‘effective temperature’ at the base 
of the extended, inflated envelope, Teff,core (see Sect. still in¬ 
creases with mass in the whole considered mass range. 


3.3. The spectroscopic HR diagram 


Figure |5] shows the location of ou r analysed models in the 
spectroscopic HR diagram (sHRD) (IFanger & Kudritzkill2014l : 
iKohler et al.ll2015h where instead of the luminosity, .if T^^lg 
is plotted as a function of the effective temperature. The quan¬ 
tity .if can be measured for stars without knowing their distance. 
Moreover, we have log(.if/jJfo) = H (cf.. Sect. 1), such that JZ” 
is directly proportional to the Eddington factor Eg as 


Fg = 


KqL 

AncGM 


eg 


Kgcr 




( 9 ) 


where g is the surface gravity and the constants have their usual 
meaning. Therefore one can directly read off Fg (right Y-axis 
in Fig. |5]) from the sHRD. Massive stellar models often evolve 
with a slowly increasing luminosity over their main-sequence 
lifetimes. Therefore, where in a conventional HR diagram mod¬ 
els with very different Fmax might cluster together (see Fig. |3, 
they separate out nicely in the sHRD since .jZ” oc LjM. The effect 
of the opacity peaks on the maximum Eddington factor (Fmax) at 
temperatures corresponding to the three partial ionization zones 
(Fe, Hell and H) is seen more clearly in the sHRD in Fig.|5]com- 
pared to the ordinary HR diagram (Fig.|2]l. 

We find that for our ZAMS models the electron scattering 
opacity is » 0.34 while the true photospheric opacity 
is around 0.5. Therefore it is expected that the true Eddington 
limit (F = 1) is achieved at about Fg = 0.7 for stellar models 
which retain the initial hydrogen abundance at the photosphere. 
Therefore in Fig. |5] we have drawn two horizontal lines corre¬ 
sponding to Fg = 0.7, one assuming the initial hydrogen mass 
fraction X - 0.74 (green line) and the other assuming Y = 0 
(red line). While models with helium-enriched photospheres ex¬ 
ceed the green line comfortably, even the most helium-enriched 
models (r otating or othe r wise) stay below the red line. 

From iKohler et al.l (1201 5h . we know that models with 
log.if/.Zf 0 > 4.4 are all hydrogen-deficient, either due to mass 
loss or due to rotationally induced mixing, as both processes 
lead to an increasing L/M-ratio (cf., their hg. 18). Figure |5]thus 
demonstrates that the models which contain super-Eddington 
layers due to the partial ionization of helium all have hydro¬ 
gen deheient envelopes, i.e., they are correspondingly helium- 
enriched. 

Figure |5] reveals that the electron scattering Eddington fac¬ 
tor Fg is not a good proxy for the maximum true Eddington 
factor (Fmax) obtained inside the star. For example, along the 
horizontal line logjZ’/.jZ^ = 4.3, corresponding to Fg ^ 0.5, 
Fmax varies from well below one to values near seven at the cool 
end. However, we note that below 30 000 K helium and hydro¬ 
gen recombine, and the gas is not fully ionised any more. The 
line opacities of helium and hydrogen become important, which 
causes the increase in Fmax (see Figs.l^andflTt. 


6 


















D. Sanyal et al.: Massive main sequence stars evolving at the Eddington limit 



Fig. 8. Maximum Eddington factor (Emax > 0.9) and surface 
Eddington factor (r(R*) > 0.9) as a function of the effective 
temperature of our analysed models. The black dotted line 
at r = 1 is drawn for convenience. 

3.4. Surface Eddington factors and the location of Emax 

The optical depth where E^ax is reached gives an idea of how 
deep in the stellar interior the layer with the highest Eddington 
factor is located. We investigate this in Fig. |6] which shows that 
Emax is located at largely different optical depths in different 
types of models. While the maximum Eddington factors oc¬ 
cur generally at optical depths below ~ 10000, we see that in 
the three effective temperature regimes identified by the super- 
Eddington peaks in Figs.|3and|3] Emax can even be located at an 
optical depth of ~ 10 or below. 

For example, when the tracks above logL/L© = 6.2 ap¬ 
proach effective temperatures of ~ 30 kK, Emax is located at the 
Fe-peak which is deep inside the envelope (t ^ 1400). The mod¬ 
els at this stage become helium-rich (7^ > 70%) and the Wolf- 
Rayet mass loss prescription is applied. Once these tracks turn 
bluewards in the HR diagram, the position of Emax jumps to the 
helium opacity peak, which is located much closer to the stellar 
surface. Consequently, we hnd three orders of magnitude of dif¬ 
ference between these two types of models with similar effective 
temperature and luminosity. 

When considering the surface Eddington-factors in the spec¬ 
troscopic Hertzsprung-Russell diagram (Fig.|2]l, we see that only 
the models with E(R*) > 0.98 have log.if/.ifj 3 values of more 
than 4.6. As discussed above, these models, which started on the 
main sequence with initial masses above 300 Mq are extremely 
helium-rich and may correspo nd to the most extreme late-type 
WNL stars (ISander et al.ll2014h . As shown in Fig.0 they exceed 
the Eddington limit by just a few per mill, which is possible be¬ 
cause of the high assumed mass loss rates that imply a slight 
deviation from hydrostatic equilibrium near the stellar surface 
(cf. Sect. [3] above). However, these models are the ones where 
our assu mption of an optica lly thin wind might break down (see 
Fig. 7 in iKohler et aLllTOT^ . Since the inclusion of an optically 
thick outflow may lead to changes of the temperature and den¬ 
sity structure near the surface, the surface Eddington-factors for 
these particular models are not reliable. 

In summary, we hnd on one hand that many of our models 
contain layers at optical depths between a few and a few thou¬ 
sand in which the Eddington factor exceeds the critical value of 
one. On the other hand, for none of our models we can con¬ 
clude that the Eddington limit is reached very near to, or at the 



Fig. 9. Density prohle of a non-rotating 285 Mq star with Eefr = 
46600 K and log(L/Lo) = 6.8 (cf. Fig.flland ApnendixlAll show¬ 
ing an inhated envelope and a density inversion. The X-axis has 
been scaled with the core radius rcore of 25.3 R©, as dehned in 
Sect. |4] 

surface, where for the vast majority we can even exclude that 
this happens. This finding leads to a shift in the expectation of 
the response of stars that reach the Eddington limit during their 
evolution. We might not expect direct outflows driven by super- 
Eddington luminosities, but instead internal structural changes, 
in particular envelope inflation. 

4. Envelope inflation 

Inflation of massive, luminous stars refers to the formation of ex¬ 
tended, extremely dilute stellar envelopes. An example of an in¬ 
flated model is shown in Fig.|3 The red shaded region is the non- 
inflated core and the blue shaded region is what we refer to as the 
inflated envelope. In the example, the model is inflated by 60% 
of its core radius (dehned below). In the presented model, the 
inhated envelope only contain a small fraction of a solar mass, 
i.e. Si 10-5 Mq. 

Envelope inhation is inherently different from classical red 
supergiant formation. The latter occurs after core hydrogen ex¬ 
haustion, as a consequence of vigorous hydrogen shell burn¬ 
ing. This process expands all layers above the shell source, 
which usually comprise several solar masses in massive stars, 
and it also operates in low mass stars, such that no proxim¬ 
ity to the Eddington limit is required. The mechanism of en¬ 
velope inhation which we discuss here works already during 
core hydrogen burning, i.e., even on the zero-age main se¬ 
quence for sufficiently luminous stars (cf., Fig.|4]l. Previous in¬ 
vestigations have suggested that inhation is related to the prox- 
imity of the ste l lar luminosity to the Eddington lumino sity 
(llshii et al.lll999t iPetrovic et al.ll2005 iGrafener et al.ll2012h in 
the envelopes of massive stars with a high luminosity-to-mass ra¬ 
tio {> 10"^Lo/Mq). The amount of mass contained in an inhated 
envelope is usually very small. As we shall see below, inhation, 
in extreme cases, can also produce core hydrogen burning red 
supergiants. 

We dehne inhation in our models through Ar/rcore : = (R*- 
rcore)/rcore, with Tcore being the radius at which inhation starts 
and Ri,, the photospheric radius. Since the densities in inhated 
envelopes are small, the dominance of radiation pressure in these 
envelopes is much larger than it is in the main stellar body. We 
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Fig. 10. Hertzprung-Russell diagram showing all the core¬ 
hydrogen burning inflated models, i.e with Ar/rcore > 0. Models 
with Ar/rcore > 5 are colored yellow. Below the solid black line 
we do not find any inflated models and above the dotted black 
line we do not find any non-inflated models in o ur grid. The 
hot pa rt of the S Dor instability strip is also marked (ISmith et al.l 
I2nn4 . 


define a model to be inflated if y6(r), which is the ratio of gas pres¬ 
sure to total pressure, reaches a value below 0.15 in the interior 
of a model. The radius at which (3 goes below 0.15 for the first 
time from the center outwards is denoted as rcore, i e- the start 
of the inflated region. The remaining extent of the star until the 
photosphere (R* - rcore) is considered as the inflated envelope. 

We emphasize that our choice of the threshold value for (3 
is arbitrary and not derivable from first principles. However, we 
have verified that this prescription identifies inflated stars in dif¬ 
ferent parts of the HR diagram very well (cf. Appendix iDl). As 
jS —> 0 for M —> oo, our criterion may fail for extreme masses. 
However, the mass averaged value of (3 for the most massive 
zero-age main sequence model analysed in the present study 
(500 Mq) is 0.3. A threshold value of 0.15 thus appears adequate 
for the present study. As an example, let us consider a typical 
inflated model, shown in AppendixlAl The value of (3 in Fig. lA.4l 
decreases sharply at the base of the inflated envelope, to around 
0.01. Even if the (3 threshold is varied by 30%, i.e. 0.15 + 0.045, 
the non-inflated core radius rcore changes by only 4%. This goes 
to show that for clearly inflated models, the value of rcore is in¬ 
sensitive to the threshold value of (3. 

We furthermore performed a numerical experiment which is 
suited to show that the core radii identified as described above 
are indeed robust. We chose an inflated 300 Mq model, and then 
increased the mixing length parameter a such that convection 
becomes more and more efficient. As shown in Fig. IBTI as a 
result the extent of the inflated envelope decreased without af¬ 
fecting the model structure inside the core radius, which thus re¬ 
mained independent of a. For a = 40 convection became nearly 
adiabatic, inflation almost disappeared, and the fact that the pho- 
tospheric radius in this case became very close to the core radius 
validates our method of identifying rcore- 

Figure [TO] shows the amount of inflation as defined above, 
for all our models that fulfil the inflation criterion, in the HR 
diagram. It reveals that overall, inflation is larger for cooler tem¬ 
peratures. This is not surprising, since inflation appears not to 
change the stellar luminosity and must therefore induce smaller 
surface temperatures. We also see that inflation is larger for more 


Fig. 11. Inflation as a function of the effective temperature for all 
analysed models which fulfill our inflation criterion. 


luminous stars, which is expected because the Eddington limit 
is supposed to play a role (see below). We also find inflation 
to increase along the evolutionary tracks of the most massive 
stars which turn back from the blue supergiant stage, which in 
this case is due to the shrinking of their core radii. A distinc¬ 
tion between the inflated and the non-inflated models is made 
by drawing the black lines in Fig. [10] They are drawn such that 
below the solid line no model is inflated and above the dotted 
line all the models are inflated. In between these two lines we 
find a mixture of both inflated and non-inflated models. We find 
that essentially all models above log(L/Lo) - 5.6 are inflated. 
Consequently, stars above ~ 40 Mq do inflate during their main 
sequence evolution. 

Figure [TT] shows the inflation factor as function of the stellar 
effective temperature for our inflated models. Whereas inflation 
increases the radius of our hot stars by up to a factor of 5, the 
cool supergiant models can be inflated by a factor of up to 40. 
We refer to Appendix for the detailed structures of several 
inflated models. 

In Fig.[T2| we take a look at inflation as a function of the core 
effective temperature Teff^core defined as 


Teff.c 


47rcrr2 


( 10 ) 


where L refers to the surface luminosity and the constants have 
their usual meaning. We can see that even our coolest models 
have high core effective temperatures, in the sense that if their in¬ 
flated envelopes were absent, their stellar effective temperatures 
would have been higher than 20 000 K. Those stars which have 
stellar effective temperatures below ~50 000 K contain the He II 
ionization zone within their envelopes, and stars with stellar ef¬ 
fective temperatures below ~10000K also contain the H/Hel 
ionization zone. However, as revealed by the density and tem¬ 
perature structure of these models (cf.. Appendix IdTi. the tem¬ 
perature at the bottom of the inflated envelope is always about 
170 000 K, and thus corresponds to the temperature of the iron 
opacity peak. We conclude that the iron opacity is at least in part 
driving the inflation of all the stars. For those with cool enough 
envelopes, helium and hydrogen are likely relevant in addition. 
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Fig. 12. Inflation as a function of the effective temperature at the 
core-envelope boundary for all analysed models with Ar/rcore > 
0. Color coding indicates the at the photosphere. 


4.1. Why do stellar envelopes Inflate? 


As suggested earlier, the physical cause of inflation in a given 
star may be its proximity to the Eddington limit. FigurefOlshows 
the correlation between inflation and Emax for our models. As 
expected, we And that our stellar models are not inflated when 
Fmax is significantly below 1, and they are all inflated for Emax > 
1. Indeed, the top panel of Fig. [T3] gives the clear message that 
the Eddington limit, in the way it is defined in Sect.[3l is likely 
connected with envelope inflation. 

Comparing Fig.[T3](top panel) to Fig.ITTIshows that inflation 
increases up to Eetr a; 5 500 K. Thereafter, and Emax decrease 
and the stars keep getting bigger without significant changes in 
Tcore, and hence, inflation still increases. However, the drop in 
inflation for the coolest models shows an opposite trend. This is 
because the non-inflated core radius rcore now moves outwards 
(increases) such that inflation (Ar/rcore) decreases even though 
/?* keeps increasing (cf. definition of rcore in Sec.lH). 

In the zoom-in at the lower panel of Fig. [T3] we see some 
models being inflated for Emax in the range ~ 0.9 ... 1. Partly, 
this may be due to the arbitrariness in our definition of inflation. 
The exact value of Ar/rcore depends somewhat on the choice of 
the threshold value of p to characterize inflation (cf. Appendix 
iDl l. i.e., the models with Arfrcom ^ 2 and Emax < 1 niay be at 
the borderline of inflation. The models with Ar/rcoie ^ 2 but 
Fmax > 1 are all very hot (Teff > 40 000 K) and in those models, 
the inflation is intrinsically small, but generally unambiguous. 

Still, we see a significant number of models below the 
Eddington limit (Emax < 1) which show a quite prominent infla¬ 
tion, i.e. which have a radius increase due to inflation of more 
than a factor of five. We investigated such a model by artifl- 
cially increasing its m ass loss rate above the critical value Merit 
(iPetrovic et al.l [2()()^ . such that the inflated envelope was re¬ 
moved (cf. Sect. l4.2b . We then found that, on turning down the 
mass loss rate to its original value, the model regained its ini¬ 
tial inflated structure with Emax < 1- However, Emax = 1 was 
reached and exceeded in the course of our experiment. We con¬ 
clude that a stellar envelope may remain inflated even if the con¬ 
dition Emax = 1 is not met any more in the course of evolution, 
but that Emax ^ 1 may be required to obtain inflation in the first 
place. 

We see that in contrast to earlier ideas of a hydrody¬ 
namic outflow being triggered when the stellar surface reaches 
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Fig. 13. Top: Inflation (Ar/rcore) as a function of Emax for the 
analysed models, the effective temperature of every model is 
color-coded. The area within the black dotted lines is magnified 
below. Bottom: Zoomed-in view of the dotted region in the top 
panel around Emax = 1 ■ 


the Eddington limit (lEddingtonl 1 1 926t lOwocki et al.l 1200^ . in 
our models this never happens, but when the properly de¬ 
fined Eddington limit is reached inside the envelope, its out¬ 
ermost layers expands hydrostatically and produce inflation. 
Two possibilities arise in this process. When the star ap¬ 
proaches the Eddington limit, the ensuing envelope expansion 
leads to changes in the temperature and the density structure. 
Consequently, the envelope opacity can either increase or de¬ 
crease. Fig.[T4]shows that the effect of expansion generally leads 
to a reduced opacity such that the expansion is indeed alleviating 
the problem. The star will then expand until the Eddington limit 
is just not exceeded any more, which is the reason why we And 
so many inflated models with Emax - 1 ■ 

Figure[T4]shows the OPAL opacities for hydrogen-rich com¬ 
position for various constant va lues of R as function of tem- 
perature, where R - p/{T /10^1^. iKippenhahn & Wei gerd (1 1 990h 
showed that for constant p - Pgas/P and constant chemical com¬ 
position, as a function of spatial co-ordinate inside the star is a 
constant. Thus, for un-inflated models, the opacity curves in Fig. 
[14] may closely represent the true run of opacity with tempera¬ 
ture inside the star. In the inflated models, p is dropping abruptly 
at the base of the inflated envelope, which means that the opac¬ 
ity is jumping from a curve with a higher /?-value to one with 
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Fig. 14. Opacity as a function of temperature for fixed values of 
the opacity parameter R defined as R = p/T^ where Fg is the 
temperature in units of 10® Kelvin. The values of \ogR which 
are held constant are indicated along the curves. The data is 
taken from the OPA L tables with a composit ion of X=0.7000, 
Y=0.2960, Z=0.004 (lldesias & Rogersl[l9^ . The black hori¬ 
zontal line shows the electron scattering opacity for a hydro¬ 
gen mass fraction of X = 0.7. At a temperature of 200 000 K, 
log/? = -5 implies a density of p = 8 x 10“*gcm“^ and 
log/? = -8 implies a density ofp = 8 x 10 " gcm“^ . 


a lower /?-value at this location. That is, the opacity is smaller 
everywhere in the inflated envelope compared to the situation 
where inflation would not have happened. 

For the chemical composition given in Fig.[T4]and assuming 
P = const., we find 


1-8 10-®y^, (11) 

such that if p drops from 0.5 in the bulk of the star to 0.1 in 
the inflated envelope, R drops by one order of magnitude. The 
corresponding reduction in opacity can be significant, i.e., up to 
about a factor of two. 

When upon expansion the envelope becomes cool enough 
for another opacity bump to come into play, the problem of not 
exceeding the Eddington limit might not be solvable this way. 
Instead, when a new opacity peak is encountered in the outer 
part of the envelope, super-Eddington conditions occur, i.e., lay¬ 
ers with Fmax > 1 (cf., Eigs.|2]and|5]), along with a strong pos¬ 
itive gas pressure (and density) gradient (cf. Sects. [2 and |5]). 
This is most extreme when the envelopes become cool enough 
(Teff < 8000 K) such that the hydrogen ionization zone is present 
in the outer part of the envelope, where Eddington factors of up 
to seven are achieved. 

4.2. Influence of mass loss on inflation 

One might wonder about the sustainability of the inflated layers 
against mass loss which is an important factor i n the evolution 
of metal-rich massive stars. iPetrovic et aP (l2006l) estimated that 
the inflated envelope can not be replenished when the mass loss 
rate exceeds a critical value of 
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Fig. 15. Mass loss history of four non-rotating evolutionary se¬ 
quences from our grid. The initial masses are given along each 
evolutionary track. The colours indicate the different mass loss 
prescriptions that were used in different phases, as described in 
Sect. 2. 


where M and Kcok stand for the stellar mass and the un¬ 
inflated radius respectiv ely, and Pniin is the m inimum density 
in the inflated region. IPetrovic et al.l (l2006ll found Merit ~ 
lO^^Moyr^* for a massive hydrogen-free Wolf-Rayet star of 
24 Mq. However, for a typical inflated massive star on the main 
sequence (see Eig. IAHT i. this critical mass loss rate is of the 
order 10“^... 10“* Moyr^'. Such high mass loss rates are ex¬ 
pected only in LBV-type giant eruptions. The mass loss rates 
app lied to our models are several orders of magnitude smaller 
(cf. lKohler etaPIMsh . 


The mass loss history of four evolutionary sequences without 
rotation are shown in Eig.[T5] We can see that even the 500 Mq 
model never exceeds a mass loss rate ~ 5 x lO^^Moyr^'. The 
critical mass loss rate for all models shown in Eig. [15] is much 
higher than the actual mass loss rates applied. Whereas Merit typ¬ 
ically exceeds M by a factor of 1000 for the inflated models in 
the 50 Mq sequence, it exceeds that of the 500 Mq sequence by a 
factor of 3 ... 100. It is thus not expected that mass loss prevents 
the formation of the inflated envelopes in massive stars near the 
Eddington limit. In fact, it may be difficult to identi fy a source of 
momentum that might drive such strong m ass loss (IShavivl200ll: 


Owocki et al.l2004 ). Grafener et al.l (1201 ih in the Milky Way and 
Bestenlehner et al.l(l2014h in the EMC found a steep dependence 


of the mass loss rates on the electron-scattering Eddington factor 
Fe for very massive stars, but they do not And mass loss rates that 
substantially exceed 10“"^Moyr“'. 


As many of the models analysed here may be pulsation- 
ally unstab l e, the mass loss rates may be enhanced in this case. 
iGrott et ^ (I 2 OO 5 I 1 show that hot stars near the Eddington limit 
may undergo mass loss due to pulsations, although extreme mass 
loss rates ar e not predicted. Eor very massive cool stars on the 
other hand, iMoriva & Langed (1201 5h And that pulsations may 
enhance the mass loss rate to values of the order of 10“^ Mq yr'^. 
Such extreme values could prevent the corresponding stars to 
spend a long time on the cool side of the Humphreys-Davidson 
limit. A detailed consideration of this issue is beyond the scope 
of the present paper. 
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Fig. 16. Upper HR diagram showing all main sequence mod¬ 
els with density inversions. Models with Ap/p > 2 have been 
colored yellow. The models without density inversion are indi¬ 
cated with open green triangles. Some representative evolution¬ 
ary tracks of non-rotating models, for different initial masses (in¬ 
dicated along the tracks in units of solar mass), are also shown 
with solid coloured lines. 


5. Density inversions 

An inflated envelope can be associated with a ‘density inversion’ 
near the stellar surface, i.e., a region where the density increases 
outward. An example is shown in Fig. |3 In hydrostatic equilib¬ 
rium, r(r) > 1 implies > 0, and thus ^ > 0. As a con¬ 
sequence, all the models which have layers in their envelopes 
exceeding the Eddington limit show density inv ersions. The cri¬ 
terion for density i nversion can be expressed as (iJoss et all 19731 
iPaxton et alj|2013h : 


Trad ^ 
Tudd 


1 -H 


/ ^Pgas \ 
WPrad/ 


(13) 


where Pgas, Prad and p stand for the gas pressure, radiation pres¬ 
sure and density respectively. A density inversion gives an in¬ 
ward force and acts as a stabilizing agent for the inflated enve¬ 
lope. We note that in the above inequality, Pgas is assumed to 
be a function of p and T only, i.e. the mean molecular weight 
p is assumed to be constant. Density inversions might also be 
present in low-mass stars like the Sun where they are caused by 
the steep increase of p around the hydrogen recombination zone 
(cf. lErgmall97lh . 

Eigure[T 6 ]identifies our core hydrogen burning models which 
contain a density inversion. The quantity Ap/p represents the 
strength of the density inversion normalized to the minimum 
density attained in the inflated zone. We can identify three peaks 
in Ap/p at Teff/kK ~ 55, 25 and 5.5 (see also Eig. [TTll. which 
coincides exactly with the three Teff-regimes in which models 
exceed the Eddington limit (cf., Fig.O. The maximum of the 
density inversions in the three zones is related to the relative 
prominence of the three opacity bumps of Fe, Hell and H re¬ 
spectively, as shown in Fig. [m 

However, an inflated model is not necessarily accompanied 
by a density inversion. This is depicted clearly in Fig. [TSl where 
we investigate the correlation between inflation and density in¬ 
version (this can also be seen by comparing Fig.[T 0 ]to Fig.fThl). 


Fig. 17. The extent of density inversions (Ap/p) as a function 
of Teff for our models. The three peaks correspond to the three 
opacity bumps of Fe, Hell and H in the OPAL tables, as indi¬ 
cated. 


Eigure [TSl shows many models which are even substantially in¬ 
flated but do not develop a density inversion. The three peaks 
in the distribution of density inversions of Fig. [17] also show up 
distinctly in this plot at the three characteristic effective temper¬ 
atures (shown in color). Models which do show a density inver¬ 
sion do always show some inflation. This is less obvious from 
Fig . [TSj because the hottest models show the smallest amount of 
inflation (Fig.fTTI). 

The stability of density inversions in stellar envelopes has 
been a matter of debate for the last few decades but th ere has 
been no consensus on this issu e yet (see Maedeij[l992h . There 
have been early speculations bv iMihal^ ( 1969h while studying 
red supergiants that a density inversion might lead to Rayleigh- 
Taylor instabilities (RTI) resulting in “elephant trunk” structures 
washing out t he positive densi ty gradient. However, as rightly 
pointed out bv ISchaereJ (Il996l) . RTI will not develop since the 
effective gravity ^eft = 0(1 - F) acting on the fluid elements is 
directed outwards in the super-Edd ington layers which contain 
the density inversion. iKutted (Il970h on the other hand claimed 
that a hydrodynamic treatment of the stellar structure equations 
will prevent any density inversion and would instead lead to 
a steady mass outflow. However, this claim is refuted by the 
present work, since our code solves the 1-D hydrodynamic stel¬ 
lar structure eq uations, in agreement with p revio us hydrodynam - 
ical models by Glatzel & Kiriakidli (Il993h and lMevne^ (Il992l) . 
IStothers & Chinll973h suggested that density inversions will 
lead to strong turbulent motions instead of drastic mass loss 
episodes. However, these layers are unstable to convection, so 
turbulence is present in any case. 

Additionally. iGlatzel & KiriakidisI (Il993h argued in favor of 
a sustainability of density inversions in the sense that they can 
be viewed as a natural consequence of strongly non-adiabatic 
convection, and they pointed out that the only plausible way to 
suppress density inversions is to use a different theory of con¬ 
vection. The only instability expected fro m simple argum ents 
therefore is co nvection which is in line with lWentzell (Il970h and 
lLangeil(ll997h . 

Still, lEkstrom et al.l (1201 2h and lYusof et ^ (1201 3h recently 
considered density inversions as ‘unphysicaT. Density inver- 
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Fig. 18. Top: Correlation between inflation and density inver¬ 
sions for our models with logL/Lo > 5. Bottom: Zoomed-in 
view of the area within the back dotted lines in the top panel. 


sions have been suppressed in their models by replacing the pres¬ 
sure scale height in the Mixing-Length Theory with the density 
scale height (cf. Sect.|7]i, as done by many stellar modelers in 
the past, often to prevent numerical difficulties. As tfie density 
scale heigfit tends to infinity when a density inversion starts to 
develop, this measure tends to enormously increase the convec¬ 
tive flux in the relevant layers. It is doubtful whether in reality 
the convective flux can be increased so much, as the ratio of the 
local thermal to the local dynamical time scale in the relevant 
layers is much smaller than one, such that convective eddies lose 
their thermal energy much faster than they rise, and thus hardly 
transport any energy at all. Multi-dimensional hydrodynamical 
simulations are desirable to settle this issue. We briefly return to 
this point in Sect. |2l 


6. Sub-surface convection 

We also studied the convective velocities in the sub-surface con¬ 
vection zones associated with the opacity peaks in our stellar 
models (ICantiello et al.ll200^ . We measure these velocities in 
units of either the isothermal or the adiabatic sound velocity, i.e. 


Cs,ad and Cs,iso respectively, which we compute as 



where y is the adiabatic index, P is the total pressure, p is the 
density, p is the mean molecular weight, T is the temperature 
and kg is the Boltzmann constant. We define Miso as the maxi¬ 
mum ratio of the convective velocity over the isothermal sound 
speed in the stellar envelope, and Mad correspondingly using the 
adiabatic sound speed. 

The true sound speed will be in between the adiabatic and 
isothermal one, closer to the first one in the inner parts of the 
star, and closer to the second in the inflated stellar envelope 
(cf.. Sect. 5). In Figs.[T9] and |20l we show the values of Miso 
and Mad for our models in the HR diagram. Whereas the convec¬ 
tive velocities are always smaller than the adiabatic sound speed, 
Figm shows that the isothermal sound speed can be exceeded 
locally in our models by a factor of a few. The convective veloc¬ 
ity and sound speed profiles for an extreme model are presented 
in Appendix ICl 

Supersonic convective velocities (adiabatic or isothermal, 
depending on the physical conditions in the envelope) may not 
be realistic and are outside the frame of the standard Mixing 
Length Theory. Therefore, in some of our models, the convec¬ 
tive velocities, and thus the convective energy transport, may 
have been overestimated. A limitation of the velocities to the 
adequate sound speed is expected to reduce the convective flux, 
which might lead to further inflation of the stellar envelope. 

The cool models which have the strongest inflation have rel¬ 
atively smaller values of Miso (compared to the hot WR-type 
models) but large values of Mad (Fig.l20l) in the sub-surface con¬ 
vection zones. This is primarily because of the fact that while 
Cs.ad depends on the total pressure, Cs,iso depends on the gas pres¬ 
sure only. In the very outer layers of the cool, luminous models, 
—> 1 and hence Pgas ~ Ptot- In such situations, Cs,ad and Cs,iso 
are only a factor sfy apart, where y is the adiabatic index. 

We And that the convective energy transport is not always 
negligible in the inflated models (cf. Sec 14. II) . We therefore eval¬ 
uate the amount of flux that is actually carried by convection 
in the inflated envelopes of our models. We define the quantity 
?7 (Mso) as the fraction of the total flux carried by convection in 
the stellar envelope, at the location where the isothermal Mach 
number is the largest. This quantity is plotted as a function of the 
effective temperature in Fig.|2T] It is evident from this figure that 
piMiso) needs not to be small for stellar envelopes to be inflated. 
However, the hotter a model is the lower its ? 7 (Miso)-value at a 
given luminosity (see Fig. l22l) . For models hotter than » 63 
kK (for e.g. the hydrogen-free He stars), Tj(Miso) indeed goes to¬ 
wards zero (Grassitelli et al., in preparation). The behaviour of 
the quantity rjiMi^o) in the HR diagram is shown in Fig. |22] 

7. Comparison with previous studies 

7.1. Stellar atmosphere and wind models 

Since the Eddington limit was thought to be reached in 
massive stars near their surface (cf.. Sect. 0, several papers 
have investigated this using stellar atmosphere calculations. 
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Fig. 19. Upper HR diagram (log L/Lq > 5) showing the maxi¬ 
mum of the ratio of the convective velocity to the local isother¬ 
mal sound speed (Miso) of the analysed models as coloured dots. 
Models with Miso < 1 have been coloured black. Some represen¬ 
tative evolutionary tracks of non-rotating models, for different 
initial masses (indicated along the tracks in units of solar mass), 
are also shown with solid black lines. 
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Fig. 21, Convective efficiency /^(Miso), which is the ratio of the 
convective flux to the total flux at the position where the isother¬ 
mal Mach number is the largest in the stellar envelope, as a func¬ 
tion of the effective temperature for all analysed stellar models 
in our grid. 
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Fig. 20. Upper HR diagram (logL/Lo > 5) showing the max¬ 
imum of the ratio of the convective velocity to the local adia¬ 
batic sound speed Mad of the analysed models (computed with 
an adiabatic index of 5/3) as coloured dots. Some representative 
evolutionary tracks of non-rotating models, for different initial 
masses (indicated along the tracks in units of solar mass), are 
also shown with solid black lines. 


Fig. 22. Upper HR diagram showing all the inflated stellar mod¬ 
els in our grid as coloured dots. The convective efficiency 
^(Miso), which is the ratio of the convective flux to the total flux 
at the position where the isothermal Mach number is the largest 
in the stellar envelope, has been color coded. Some represen¬ 
tative evolutionary tracks of non-rotating models, for different 
initial masses (indicated along the tracks in units of solar mass), 
are also shown with solid black lines. 


iLamers & Fitzoatrickl (Il988l) investigated the Eddington fac¬ 
tors in the atmospheres of luminous stars in the te mperature - 
gravity diagram, while lUlmer & Fitzpatri^ (Il998h did so in 
the HR diagram. Both studies took the full radiative opacity 
into account. While their technique did not allow them to reach 
Eddington factors of one or more, Ulmer and Fitzpatrick found 
that model atmospheres with a maximum Eddington factor of 
0.9 are located near the observe d upper luminosity limit of stars 
in the Large Magellanic Clou d (iHumphrevs & Davidsonl[r979t 
iFitzpatrick & Garmanvll9^ . 


The main feature in the lines of constant Eddington-factors 
in the HR diagram found by Ulmer and Fitzpatrick is a drop 
from Teff ^ 60 000 K to 15 000 K. This may correspond to the 
drop in the maximum Eddington factor seen in our models in the 
same temperature interval (cf., Figs.|2]and|2l- Note that the peak 
around Tefr ^ 30 000 K in Fig.[3]corresponds only to helium-rich 
models, which are not considered by Ulmer and Fitzpatrick. 

On the other hand, neither inflation nor super-Eddington lay¬ 
ers or density inversions are reported by Ulmer and Fitzpatrick, 
or from any hot, main sequence star model atmosphere calcula¬ 
tion so far (to the best of our knowledge). One reason might be 
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that many model atmospheres only include a rather limited op¬ 
tical depth range (e.g., up to t = 100 in Ulmer and Fitzpatrick), 
such that the iron opacity peak is often not included in the model. 
Additionally, the computational methods employed might not al¬ 
low for a non-monotonic density profile. 

Given the ubiquity of inflation for models above log L/Lq > 
5.5 or M > 50 Mq in the LMC, and a correspondingly lower 
limit in the Milky Way due to its higher iron content, it is de¬ 
sirable to construct model atmospheres which include this effect 
and identify its observational signatures. As the density profiles 
of such atmospheres near the photosphere are significantly dif¬ 
ferent from those in non-inflated atmospheres, such signatures 
ma y indeed be expe cted. 

lAsnlundl (Il998h gives a thorough analysis of the Eddington 
limit in cool star atmosphere models. He does indeed find super- 
Eddington layers and density inversions in his models, and gives 
arguments for the physically appropriate nature of these phe¬ 
nomena. He also discusses the effects of stellar winds on these 
features, and finds they may be suppressed by extremely strong 
winds, but not by winds with mass loss rates in the observed 
range. Asplund does not find inflation in his models, arguably 
because again, the iron opacity peak is not included in his model 
atmospheres, which appears essential even for our models with 
co ol effective temperatu res. 

lOwocki et al.l il2004ll and Ivan Marie et all (l2008l) studied the 
winds of stars which reach or exceed the Eddington limit at 
their surface. As we have shown above, this condition is gen¬ 
erally not found in our models (cf., Eigs.|7]and|3l. However, it 
may occur in helium-rich s tars (see again Eig.|7]l an d hydrogen- 
free Wolf-Rayet stars (cf. iHeeer & Langer 1996h . as well as 
in stars whi ch deviate from ther mal or hydrostatic equilibrium. 
Noticeably, lOwocki et al.l (l2004l) find that the mass loss rates in 
this case are still quite limited, due to the energy loss attributed 
to lifting the wind mate rial out of the gravitational potential (see, 
iHeger & Langeilfl996l) . 

We want to emphasize in this context that the Eddington 
limit investigated in the quoted models as well as in our 
own may be different from the true Eddington limit, due to 
a number of effects which are all related to the opacity of 
the stellar matter in the stellar envelopes. One is that con¬ 
vection, which is necessarily present in the layers near or 
above the Eddington limit, may induce density inhomogeneities 
or clumping which can alter the effective radiative opacity 
(HE avivlll998h . In fact, depending on the nature of the clu mp- 
ing, the opacity may be en h anced (Grafener et al.l l2012h or 
reduced jOwocki et al.l 120041: iRuszkowski & BegelmanI 1^031 
iMuiires et al. 20111) . Eurthermore. such opacity calculations are 
tedious, and even in the currently used opacities, important con¬ 
tributions might still be missing. 

Einally, the effect of stellar rotation on the stability limit 
in the atmospheres especially of hot stars is clearly important 
(lLangerlll997lll99^lMaeder & Mevn^l2000l) . However, it adds 
another dimension to this difficult problem and is therefore gen¬ 
erally not included (cf. Sect. O- 


7.2. Stellar Interior models 


The peculiar core-h alo density structure of i nflated stars has first 
been pointed out bv IStothers & ChinI d 19931) . after the large iron 
bump in the opacities near 170 000 K was found bv llglesias et al.l 
19921). Eurther s t udies pointing out this phenomenon comprise 


ilshii et all ( 1999). Retro vie et al.l (l2006l) . iGrafener et al.l (1201 2h 


and 


Kohler et al 


( |2015|) . Conceivably, inflation may be present 
in further models of very massive stars, but often no statements 


on the presence or absence of this phenomenon are made in the 
respective papers. 

Eor ex a mple, the models for very massive stars by 
lYusof et aH (1201 3l) only discuss the electron-scattering 
Eddington factor in their models. On the ZAMS, their models 
are ho tter and more compact than the ones of iKohler et al.l 
(I 2 OI 5 I) analysed here, which implies that inflation is either 
weaker or absent. This difference might be due to the different 
treatm ent of conve c tion i n the sub-surface convective zones, 
where lYusof et aH (1201 3l) assume the mixing length to be 
proportional to the density scale height instead of the standard 
pressure scale height. This prohibits the formation of density 
inversions, and since the density scale height tends to infinity 
when a model attempts to establish a density inversion, the 
convection may transport an arbitrarily large energy flux in this 
scheme. While the physics of convection introduces one of the 
biggest uncertainties in the atmospheres of stars close to the 
Eddington limit, efficient convective energy transport in inflated 
envelopes appears unlikely (cf. Sect.|5]). 

A suppression of inflation may have significant conse¬ 
quences for the evolution of massive stars, as the stellar mod¬ 
els stay bluer and as a result have lower mass loss rates and 
lower spin-down rates. The final fates of such non-inflated stars 
will be significan tly different compared to inflated stars (see 
iKohler et al.ll2015L for a detailed discussion). 


Grafener et al.l (1201 2l) find inflation which, for their mod¬ 


els wit hout clumping, correspond well to those of iPetrovic et al.l 
(I 2 OO 6 I) for the Wolf-Rayet case, and to our unpublished solar 
metallicity main sequence models, which show a bending of the 
zero-age main sequence to coo l temperatures for M > 100 Mq. 
The models ofllshiLeLalJ (ll99 ^ aga in agree very well. Including 
the work of IStothers & Chini 'i 19931) . we conclude that the effect 
of inflation in models of massive main sequence stars is found 
in at least four independent stellar structure codes, with three of 
them quantitatively producing very similar results. 

As pointed out above, massive star evolutionary models 
which include effects of rotation are being produced rou¬ 
tinely these days (cf., Maeder & Mevned 1201 Ol: iLangerl 1201 2t 
IChieffi & Limongill2013l) . but an investigation of the effect of 
stellar rotation on the stability limit in the atmospheres of hot 
stars requires the construction of two-dimensional stellar mod¬ 
els. 


8. Comparison with observations 

8.1. The VFTS sample 

A prime motivation of IKohler et al.l (l2015l) for computing the 
evolutionary models for the very massive stars analysed here was 
to provide a the oretical framework for the VET Elames-Tarantula 
Survey (VETS. [Evans et aTll201 ll) . Within VETS, multi-epoch 
spectral data of about 700 early B and 300 O stars are being anal¬ 
ysed thro ugh detailed model atmos pher e calculations. Within 
this effort, iBestenlehner et al] (l2014l) and iMcEvov et al.l (l2015h 
derived the physical properties of more than 50 very massive 
stars, with luminos i ties lo g(L/ Lq) >5.5. We confront the mod¬ 
els of iKohler et al.l (l2015l) with this sample in Eig. |23] 

Two sets of model data are included in Eig. one which 
uses the effective temperatures of the Kohler et al. models di¬ 
rectly, and a second one where the effective core temperature is 
used as defined in Sect.|4](cf. Ea.fTOll. The latter approximates 
the surface temperature of our models if inflation was com¬ 
pletely absent. An example calculation presented in Appendix 
m where inflation in a 300 M 0 is suppressed by increasing the 
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mixing length parameter, shows that this approximation is in¬ 
deed quite good. The zero-age main sequence is also drawn for 
both sets of models. Note that while wind effects are clearly seen 
in the spectra of all stars in the sample, the optical depth of their 
winds is expected not to exceed t = 2 (cf., fig. 7 in iKohler et al.l 
(l2015h ) until the stars become very helium rich at their surface. 
Therefore, the effective temperatures derived from the observa¬ 
tions need not be corrected for optic a lly thi ck winds. 

As shown in iBestenlehner et al] (l2014l) . the hottest stars in 
their sample follow the ZAMS of the Kohler et al. models very 
closely, well into the regime of inflation. One might expect un¬ 
evolved stars to the left of the Kohler et al. ZAMS if inflation 
was not present. In that case, the stars above logL/Lo ^ 6.2 
might spend a significant fraction of their life time on the hot 
side of the Kohler et al. ZAMS. The absence of such hot stars, 
however, does not conclusively argue that inflation does exist in 
nature, i.e., even without inflation, the star formation history in 
30 Doradus might preclude the existence of such stars, or they 
may be hidden in their natal cloud due to their youth jYorkel 
[I^ICastro 61^12014 . 

On the cool side, Fig.l23]shows an absence of observed stars 
for logL/Lo ^ 6.15 and Tetr < 35 000 K. As the evolutionary 
models predict about 30% of the core hydrogen burning to take 
place at Tefr < 35 000 K in this luminosity regime, this may indi¬ 
cate that the inflation in the models of Kohler et al. is too strong. 
On the other hand, again, the absence of correspondingly cool 
stars 30 Doradus may be a result of the local star formation his¬ 
tory. 

In the luminosity range below, at 5.5 < logL/Lo < 
6.15, stars as cool a s Tefr < 15 000K are observed for which 
IMcEvov et al.l (l2015h concluded that they are still core hydrogen 
burning objects. The observed stars are somewhat cooler than 
the coolest core effective temperature of our models, which may 
argue in favor of inflation in real stars. Note that the life time of 
stars in the regime < 20 000 K is only 10% for the Kohler et 
al. models in the considered luminosity range. 

In summary, as the stellar evolution models for these high 
masses are still quite uncertain, we And it not possible to argue 
for or against inflation being present in the observed stars con¬ 
sidered here based on Fig.|23] In fact, it is intriguing that most of 
the observed stars are found in the regime where the inflated and 
non-inflated models overlap. Nevertheless, the observed sample 
above logL/Lo - 5.5 might constitute the best test case, since 
according to our models, the envelopes of all of them are ex¬ 
pected to be strongly affected by the Eddington limit. Model 
atmosphere calculations for these stars which include inflation 
might shed new light on this question. 

8.2. Further possible consequences of inflation 
S Doradus type variability 

iGrafener et al.l (1201 2h argue that the S Doradus type variabil¬ 
ity of LB Vs may be related to the effect of inflation, and fo - 
cussed in particular on the case of AG Car (iGroh et al.ll2009h . 
They propose that an instability sets in when their 70 Mq chem¬ 
ically homogeneous hydrostatic stellar model is highly inflated 
(=s 120 Rq) by virtue of which the inflated layer becomes gravi¬ 
tationally unbound and a mass loss episode follows. 

In contrast to this idea, our inflated hydrodynamic stellar 
models do not show any signs of such an instability. This could 
be due to various simp lifying assumptions made in the models of 
IGrafener et akl (120121) . in particular their neglect of the convec¬ 
tive flux in the inflated envelopes (cf., Sectl6]l. Nevertheless, the 


physics of convection is very complex in these envelopes, and 
our results do not imply that instabilities may not occur. 

In fact, consideri ng the hot edge of the S Doradus variabil¬ 
ity strip according to ISmith et^ (|2004 in Fig. |23 we see that 
it roughly separates the models with a low maximum convec¬ 
tive efficiency (//max < 0.2) from those with a higher convective 
efficiency. If such high fluxes would not be achievable in these 
envelopes (cf.. Sect. |6]l, a dynamical instability might well be 
possible. 

The hot edge of the S Doradus variability in the HR dia¬ 
gram also coincides quite well with the borderline separating 
mildly (Ar/rcore < 1) from strongly (Ar/rcore >1) inflated mod¬ 
els. Comparing this with the observed distribution of very mas¬ 
sive stars in Fig.|^ which indicates that essentially no stars are 
found far to the cool side of this line, could indicate again that 
strongly inflated envelopes are indeed unstable, and might lead 
to S Doradus type variability and an increased time-averaged 
mass loss rate. 


LBV eruptions 

iGlatzel & KiriakidisI (1 19931) speculated that strange mode pulsa¬ 
tions might be responsible for the LBV phenomenon. These pul¬ 
sations are characterized by very short growth times (~ Tdyn) and 
small brightness fluctuations roughly of the order ~ 10... 100 
mmag dGlatzel et alJl99^lGrott et al.l2005|) . However, the mass 
contained in the pulsating envelopes of their models is negligi¬ 
ble compared to the stellar mass, and the associated brightness 
variations cannot explain the humongous luminosity variations 
observed in LBV eruptions. 

We have seen in Sect. |5] that an inflated envelope often pro¬ 
duces a density inversion. Such density inversions have been re¬ 
peatedly proposed as a so urce of instabi l ities giving rise to erup- 
tive mass loss in LBVs (lMaedeij|1989l: iMaeder & Contiiri994 
IStothers & ChinI [19931) . Given our results, we consider it un¬ 
likely that a density inversion can be the sole cause of LBV erup¬ 
tions. Density inversions are a generic feature present in a mul¬ 
titude of our models (see Fig. [S, while the LBV phenomenon 
is quite rare. Furthermore, the density inversions in our models 
are found very close to the surface of the star, with very small 
amounts of mass above it. 

Given our results, inflation per se appears unlikely to cause 
LBV eruptions, again, because it occurs too abundantly in our 
models, and also because the mass of the inflated envelope is 
generally very small. However, Fig.|24]reveals that this is not so 
for our coolest models. Whereas for most models the mass of 
the inflated envelope is smaller than ~ 10“^ M©, intriguingly it 
rises to several solar masses in the models which have effective 
temperatures below ~ 10 000 K. These cool models, of which 
detailed examples are presented in Fig.|Dl also show the highest 
Eddington factors (Fig.[3]l and the strongest inflation (Fig.fTTb. 
This behaviour is seen in the mass range of 40... 100 M©, which 
corresponds well to the masses of observed LBVs. 

A key feature in our cool models with massive inflated en¬ 
velopes is visible in Fig.|2T] As the opacity in the hydrogen re¬ 
combination zone becomes very large, effectively blocking any 
radiation transport, convective efficiencies of the order of t; ^ 1 
are needed to transport the stellar luminosity through this zone. 
Also, in the iron convection zone at the bottom of the envelope, a 
high convective efficiency (rj ^ 0.5) is found in these models. As 
shown in Figs.[T9]and|20l this requires sonic or even supersonic 
convective velocities in the framework of the standard MET as 
implemented in our code. It thus appears conceivable that in re¬ 
ality convection is less efficient in such a situation, for e.g., be- 
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Fig. 23. Hertzsprung-Russell diagram showing the Koehler et al. models which include inflation (red) and corresponding non- 
inflated st ellar models (blue) above l og L/L© > 5 (see text for explanation). The triangles refer to all the Of/WN and W Nh single stars 
studied bv iBestenlehner et al.l(l2014h whereas the black dots and the diamonds refer to the B-supergiants (single) from lMcEvov et all 
(1201 5h and WN5h sta rs of the core R136 fro m Crowther et al.l (l2010l) . respectively. The O stars observed within the VFTS survey are 
marked with crosses (ISabm-Saniulian et al.ll2014l Ramirez-Agudelo et al. in preparation). The zero-age main-sequence of the non¬ 
rotating stars is marked with the solid line while the dashed line indicates the approximate position of t he zero-age main se quence 
if inflation was absent (see text for further details). The hot part of the S Doradus instability strip from ISmith et akl (l2004l) is also 
shown for reference. 


cause of viscous dissipation. When a star enters this region with 
a massive inflated envelope, throughout which the stellar lumi¬ 
nosity can neither be transported by radiation nor by convection, 
hydrostatic equilibrium will not be possible any more, and the 
loosely bound inflated envelope may by dynamically ejected. We 
believe that this scenario may relate to major LBV eruptions. 



Fig. 24. Mass contained in the inflated envelope for all inflated 
models of our grid, as a function of the effective temperature. 
The effective actual mass of the models is colour-coded (see 
colour bar to the right). 
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As proposed by iLangeil (1201 2h . a rapid evolutionary time 
scale of a star may be required to obtain LBV outbursts in ad¬ 
dition to the star reaching the Eddington limit. If evolution on 
the thermal timescale is rapid enough, this might produce LB Vs 
after core hydrogen exhaustion which may relate to most of the 
observed LB Vs in our Galaxy and the Magellanic Clouds, as 
well as LB Vs after core helium exhaustion, which may concern 
to the recently accumulated e vidence of LBV outbu rsts in imme¬ 
diate supernova progenitors (ISmith & Arnettll2014h . 

Supernova shock break-out 

A recent study by iMoriva et ahl (1201 5h concluded that inflated 
stellar models can help to explain the extended rise time of the 
shock break-out signa l from the Type Ib supernova SN2008D 
dSoderberg et al.ll2008l) . In this scenario the shock break-out oc¬ 
curs deep inside the inflated envelope and consequently the rise 
time is determined by the radiative diffusion time of the enve¬ 
lope and not the light crossing time. They also noted that more 
such events, if observed in future, might serve as indicators of 
inflated supernova progenitors. 

Whereas the above result supports the idea that hydrogen- 
free Wolf-Ra y et sta rs may posse ss infl ated envelopes 
(iPetrovic et al.l 120061: iGrafener et al.l 1201 2h . LB Vs have 
also been suggested to be immediate pr ogenitors of supernovae 
dKotak & Vinld [200^ iGroh et ^ l2013h . In the realm of high 
cadence supernova surveys, this opens up the possibility to also 
test the existence of envelope inflation in hydrogen-rich stars 
through supernova shock break-out observations. 

9. Discussion and conciusions 

We investigated the internal structu res of the massive star m od- 
els computed by iBrott et ^ d201 ih and iKohler et al.l d2015h us¬ 
ing a 1-D hydrodynamical stellar evolution code, with particular 
emphasis to the Eddington limit. We And that the conventional 
idea of sufficiently massive stars reaching the Eddington limit at 
the stellar surface is not reproduced by our core hydrogen burn¬ 
ing models, not even at 500 Mq (cf., Eigs.|7]and[8]l. Instead, we 
find a suitably defined Eddington limit inside the star (Eqn.[3]l is 
reached by models with log(L/ L©) >5.6 (Eig. la , which leads to 
sub-surface convection, envelope inflation (Eig.fOl) and possibly 
to pulsations. Many of our models even exceed this Eddington 
limit, in the extreme case of red supergiants even by factors of 
up to seven (Eig. [3]), with the consequence that strong density in¬ 
versions develop such that hydrostatic equilibrium is maintained 
(Eig.[ia. 

In the analysed models, whose initial composition is chosen 
to match that of the LMC, all stars above ~ 40 M© do reach the 
Eddington limit in their envelopes. As iron opacities are mainly 
responsible for this phenomenon, we expect that this mass limit 
is higher at a lower metallicity, and similarly lower for massive 
stars in our galaxy. Eurthermore, there may be two groups of 
stars for which this limit comes down even further. Eirstly, the 
centrifugal force in rapidly rotating stars may lead to similar 
conditions in the envelope layers near the stellar equator, i.e. to a 
strong latitude dependence of inflation. Perhaps, this could give 
rise to the so called B[e] supergiants, which show a slow and 
dense equatorial win d and a fast polar wind at the same time 
(IZickgraf et al.llT98^ . Second, the mass losing stars in interact¬ 
ing close binary systems ev olve to much higher L/M-va lues than 
corresponding single stars dLanger & Kudritzldl2014li . and are 
therefore expected to reach the Eddington limit for much lower 
initial masses. 


The stability of the inflated envelopes is not investigated 
here , but many of them are l i kely to be puls ationally unsta¬ 
ble (iGlatzel & Kiriakidiill993l:ISaio et al.lll998l Sanyal et al. in 
preparation). If so, it is expec ted that the pulsations w ill lead to 
mass loss enhancements (e.g. lMoriva & Langerll2015l) , or to the 
loss of the inflated envelope. In the latter case, the envelope is 
expected to re-grow unless the achieved time average mass loss 
rate exceeds the high critical mass loss rate (Sect. lOi . We find 
that in our coolest models, the mass contained in the inflated en¬ 
velopes can reach several solar masses (Fis.UM. and speculate 
that their dynamical loss may resemble LBV major eruptions (cf. 
Sect. 18.21 1. Consequently, even though reaching or exceeding the 
Eddington limit may not immediately lead to strong outflows 
in stars, clearly the mass loss rate of the stars will be strongly 
affected, in the sense that the mass loss will be significantly en¬ 
hanced one way or another. 

It will be crucial to test observationally whether luminous, 
main sequence stars indeed possess inflated stellar envelopes. 
This possibility has not yet been investigated with stellar atmo¬ 
sphere models for h ot stars. Perhaps the b est candidates are the 
S Doradus variables (IGrafener et al.ll20T^ . which appear in the 
part of the HR diagram where our models predict a radius infla¬ 
tion by more than a factor of two (cf.. Pig. [loll. 

Pinally, we note that besides massive stars, the Eddington 
limit is relevant to various other types of stars, as luminous post- 
AGB star. X-ray bursts. Novae, R Corona Borealis stars and ac¬ 
creting compact objects. It may be interesting to assess to what 
extent similar phenomena as found in this work might play a role 
in these objects. 
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Fig. A.l. Density structure of the stellar model. The black solid 
line marks the base of the inflated envelope, i.e. where /? = 0.15. 
The intersection of the dotted lines with the red line on either 
side mark the points where/? = 0.15 + 0.045. 



Fig. A.2. Run of T in the interior of the stellar model. 

Appendix A: Interior structure of a 85 Mq stellar 
model 

Appendix B: Effect of efficient convection on 
inflation 

Knowing that convective flux is proportional to the mixing 
length, we show here (Fig. IbtT i that by increasing the mixing 
length parameter a in an inflated 300 Mq model near the ZAMS, 
the inflation gradually goes away and what we are left with is an 
almost non-inflated star, whose radius is well-approximated by 
core radius rcore of the inflated model. 


Fig.A.3. Rosseland mean opacity k in the interior of the stellar 
model. 



Fig.A.4. Run of yS(= Rgas/Rtot) in the interior of the stellar 
model. 

In such models, turbulent pressure becomes important (which is 
not taken into account in our models) as well as standard MLT 
fails to be a good approximation for modelling convection. 

Appendix D: Representative models 

The profiles of different relevant physical quantities are shown 
for a few selected stellar models at five distinct effective tem¬ 
peratures corresponding to the three peaks in Tmax and the two 
troughs in between the peaks (cf., Fig.|2]). 


Appendix C: Convective velocity profile in a WR 
model 

The convective velocity is shown as a function of radius in 
a massive (147 Mq) WR-type (Tj = 0.89) stellar model, in 
Fig. 1C. II The variation of isothermal and adiabatic sound speeds 
are also plotted for comparison. The convective velocities ex¬ 
ceed the local isothermal sound speed in the envelope where con¬ 
ditions are non-adiabatic, i.e. thermal adjustment time is short. 
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Fig.A.5. Fraction of flux carried by radiation (Liad/Ltot) in the 
interior of the stellar model. 



Fig. B.l. Density profile of a 300 M© model with different val¬ 
ues of the mixing length parameter a (see Sec.|2]i. The black 
dotted line marks the location of rcore, i e. the base of the inflated 
envelope where = 0.15. 



Radius [Rq] 

Fig.C.l. Convective velocity, isothermal sound speed and adi¬ 
abatic sound speed profiles in a 147 M© WNL type star with 
Y, = 89%. 
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Fig.D.l. Detailed structure examples for stellar models with an effective temperature near 50 000 K, for three different luminosities 
(cf., Fig. 2). The dashed line marks the point at which (3 falls below 0.15, i.e. the beginning of the inflated envelope. The square 
symbol marks the temperature Tpe at which k is maximum due to the iron opacity bump. The hatched regions show the convective 
zones. 
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Fig.D.2. Detailed structure examples for stellar models with an effective temperature near 25 000 K, for three different luminosities 
(cf., Fig. 2). The dashed line marks the point at which (3 falls below 0.15, i.e. the beginning of the inflated envelope. The square and 
the cross mark the temperatures Tpe and Tpe at which k is maximum due to the iron and the helium opacity bumps respectively. The 
hatched regions show the convective zones. 
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Fig.D.3. Detailed structure examples for stellar models with an effective temperature near 5 000 K, for two different luminosities 
(cf., Fig. 2). The dashed line marks the point at which (3 falls below 0.15, i.e. the beginning of the inflated envelope. The square, 
cross and the circle mark the temperatures Tpe, Tpe and Th at which k is maximum due to the iron, helium and hydrogen opacity 
bumps respectively. The hatched regions show the convective zones. 
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Fig. D.4. Detailed structure examples for stellar models with an effective temperature near 32 000 K, for three different luminosities 
(cf., Fig. 2). The dashed line marks the point at which (3 falls below 0.15, i.e. the beginning of the inflated envelope. The square 
symbol marks the temperature Tpe at which k is maximum due to the iron opacity bump. The hatched regions show the convective 
zones. 
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Fig.D.5. Detailed structure examples for stellar models with an effective temperature near lOOOOK, for three different luminosities 
(cf., Fig. 2). The dashed line marks the point at which (3 falls below 0.15, i.e. the beginning of the inflated envelope. The square, 
cross and the circle mark the temperatures Tpe, Tpe and Th at which k is maximum due to the iron, helium and hydrogen opacity 
bumps respectively. The hatched regions show the convective zones. 
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